Code_deliniation

Author

Joël Gaxotte

Watershed Delineation

Ce document présente un workflow complet pour la délimitation de bassins versants à partir d’un Modèle Numérique d’Élévation (MNE). Le processus inclut le prétraitement du MNE, la génération de réseaux hydrographiques et la délimitation finale des bassins versants.

Chargement des packages

#Packages de base pour la manipulation de données spatiales
options(repos = c(CRAN = "https://cloud.r-project.org"))
install.packages("Rcpp")

The downloaded binary packages are in
    /var/folders/qb/stf0vsnn0zn0509g0d3ch4r40000gn/T//Rtmp04WMM7/downloaded_packages
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(raster)
Loading required package: sp

Attaching package: 'raster'

The following object is masked from 'package:dplyr':

    select
library(terra)      # Alternative moderne au package raster
terra 1.8.21

Attaching package: 'terra'

The following object is masked from 'package:tidyr':

    extract
library(sf)         # Pour les objets vectoriels
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(tmap)       # Cartographie interactive
library(stars)      # Autre format de données spatiales
Loading required package: abind
library(whitebox)   # Interface avec WhiteboxTools

# Initialisation de WhiteboxTools
whitebox::wbt_init()
wbt_version()
[1] "WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020"                  
[2] ""                                                                           
[3] "WhiteboxTools is an advanced geospatial data analysis platform developed at"
[4] "the University of Guelph's Geomorphometry and Hydrogeomatics Research "     
[5] "Group (GHRG). See https://jblindsay.github.io/ghrg/WhiteboxTools/index.html"
[6] "for more details."                                                          
# Configuration des thèmes et modes d'affichage
theme_set(theme_classic())
tmap_mode("view")
ℹ tmap mode set to "view".

Importation et préparation du MNE

Cette section importe le Modèle Numérique d’Élévation (MNE) et effectue les premiers traitements pour éliminer les valeurs aberrantes.

# Lecture du fichier MNE (format .adf)
MNE <- rast("/Users/thejo/Desktop/TD/Donnees_Audrey/DEM/demwatershed")

# Nettoyage des valeurs négatives (souvent des artefacts)
MNE[MNE < 0] <- NA

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          
# Visualisation interactive du MNE
tm_shape(MNE) +
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE) +
  tm_scale_bar()
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
`tm_scale_continuous()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values') to
  'tm_scale_continuous(<HERE>)'
! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.
SpatRaster object downsampled to 3740 by 2675 cells.

[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "PuOr" is named
"brewer.pu_or"

2. Conversion et ombrage du MNE

Conversion au format .tif pour une meilleure compatibilité et génération d’un ombrage pour la visualisation.

# Conversion en format .tif
writeRaster(MNE, "/Users/thejo/Desktop/TD/trois/donnees/MNE.tif", overwrite = TRUE)

|---------|---------|---------|---------|
=========================================
                                          
# Génération d'un ombrage (hillshade)
wbt_hillshade(
  dem = "/Users/thejo/Desktop/TD/trois/donnees/MNE.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/hillshademne.tif",
  azimuth = 115  # Angle d'éclairage en degrés
)
[1] "hillshade - Elapsed Time (excluding I/O): 12.842s"
# Visualisation de l'ombrage
hillshade2 <- rast("/Users/thejo/Desktop/TD/trois/donnees/hillshademne.tif")
tm_shape(hillshade2) +
  tm_raster(style = "cont", palette = "-Greys", legend.show = FALSE) +
  tm_scale_bar()
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
`tm_scale_continuous()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values') to
  'tm_scale_continuous(<HERE>)'
[v3->v4] `tm_raster()`: use `col.legend = tm_legend_hide()` instead of
`legend.show = FALSE`.
! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.
SpatRaster object downsampled to 3740 by 2675 cells.

[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "-Greys" is named
"greys" (in long format "brewer.greys")

Correction des dépressions

Les dépressions dans le MNE peuvent affecter l’analyse hydrologique. Cette section les corrige.

# Conversion en format float32 avec valeur NA spécifique
terra::writeRaster(
  MNE, 
  "/Users/thejo/Desktop/TD/trois/donnees/MNE_float32.tif", 
  datatype = "FLT4S", 
  overwrite = TRUE, 
  NAflag = -9999
)

|---------|---------|---------|---------|
=========================================
                                          
# Correction des dépressions par la méthode du moindre coût
wbt_breach_depressions_least_cost(
  dem = "/Users/thejo/Desktop/TD/trois/donnees/MNE_float32.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/MNE_breach_dep.tif",
  dist = 100,  # Distance maximale de correction, à ajuster selon le cas
  fill = TRUE
)
[1] "breach_depressions_least_cost - Elapsed Time (excluding I/O): 34min 54.999s"
# Remplissage complémentaire des dépressions
wbt_fill_depressions_wang_and_liu(
  dem = "/Users/thejo/Desktop/TD/trois/donnees/MNE_breach_dep.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif"
)
[1] "fill_depressions_wang_and_liu - Elapsed Time (excluding I/O): 1min 22.107s"
MNE_clean <- rast("/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif")

# Vérification des résultats
summary(values(MNE_clean))
  MNE_cleaned       
 Min.   :  0        
 1st Qu.:475        
 Median :541        
 Mean   :500        
 3rd Qu.:591        
 Max.   :944        
 NA's   :131419502  

Analyse du réseau hydrographique

Calcul de l’accumulation de flux et de la direction d’écoulement.

# Accumulation de flux (méthode D8)
wbt_d8_flow_accumulation(
  input = "/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/flow_accD8.tif"
)
[1] "d8_flow_accumulation - Elapsed Time (excluding I/O): 11.71s"
  # Transformation logarithmique pour meilleure visualisation sur ArcGIS pro
flow_acc <- rast("/Users/thejo/Desktop/TD/trois/donnees/flow_accD8.tif")
scaled_acc <- log1p(flow_acc)

|---------|---------|---------|---------|
=========================================
                                          
terra::writeRaster(
  scaled_acc, 
  "/Users/thejo/Desktop/TD/trois/donnees/flow_accD8_scaled.tif", 
  overwrite = TRUE
)

|---------|---------|---------|---------|
=========================================
                                          
# Direction d'écoulement (méthode D8)
wbt_d8_pointer(
  dem = "/Users/thejo/Desktop/TD/trois/donnees/MNE_cleaned.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/pointer3.tif"
)
[1] "d8_pointer - Elapsed Time (excluding I/O): 3.470s"
pointer3 <- rast("/Users/thejo/Desktop/TD/trois/donnees/pointer3.tif")

# Vérifier les valeurs du flow direction (Ne devrait que contenir: 1,2,4,8,16,32,64,128)
freq(pointer3)
  layer value   count
1     1     0    2605
2     1     1 4809562
3     1     2 5035026
4     1     4 3683410
5     1     8 5557440
6     1    16 4636193
7     1    32 4522006
8     1    64 3083810
9     1   128 4383217

Définition des points d’exutoire

Avec point(s) déjà défini(s)

Création et préparation des points d’exutoire pour la délimitation des bassins versants.

# Création des points d'exutoire
ppoints <- tribble(
  ~Lon, ~Lat,
  -63.369243, 50.828091
)

# Conversion en objet spatial
ppointsSP <- SpatialPoints(ppoints, proj4string = CRS("+proj=longlat +datum=WGS84"))

# Export en shapefile
shapefile(ppointsSP, filename = "pourpoints.shp", overwrite = TRUE)

Extraction du réseau hydrographique

Seuillage de l’accumulation de flux pour identifier les cours d’eau.

wbt_extract_streams(
  flow_accum = "/Users/thejo/Desktop/TD/trois/donnees/flow_accD8.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/raster_streams.tif",
  threshold = 200  # Ajuster selon la résolution et la zone d'étude
)
[1] "extract_streams - Elapsed Time (excluding I/O): 2.392s"
wbt_jenson_snap_pour_points(
  pour_pts = "/Users/thejo/Desktop/TD/pourpoints.shp",
  streams = "/Users/thejo/Desktop/TD/trois/donnees/raster_streams.tif",
  output = "/Users/thejo/Desktop/TD/trois/donnees/snappedpp.shp",
  snap_dist = 200  # Distance d'ajustement (unités du SCR)
)
[1] "jenson_snap_pour_points - Elapsed Time (excluding I/O): 0.2s"
snappp <- shapefile("/Users/thejo/Desktop/TD/trois/donnees/snappedpp.shp")
streams <- raster("/Users/thejo/Desktop/TD/trois/donnees/raster_streams.tif")

#Vérifier concordance des systèmes de coordonnées
print("Pointer3 CRS:")
[1] "Pointer3 CRS:"
print(crs(pointer3))
[1] "PROJCRS[\"unnamed\",\n    BASEGEOGCRS[\"Unknown datum based upon the GRS 1980 ellipsoid\",\n        DATUM[\"Not specified (based on GRS 1980 ellipsoid)\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101004,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4019]],\n    CONVERSION[\"Lambert Conic Conformal (2SP)\",\n        METHOD[\"Lambert Conic Conformal (2SP)\",\n            ID[\"EPSG\",9802]],\n        PARAMETER[\"Latitude of false origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8821]],\n        PARAMETER[\"Longitude of false origin\",-95,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8822]],\n        PARAMETER[\"Latitude of 1st standard parallel\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8823]],\n        PARAMETER[\"Latitude of 2nd standard parallel\",77,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8824]],\n        PARAMETER[\"Easting at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8826]],\n        PARAMETER[\"Northing at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"easting\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"northing\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"
print("Snappp CRS:")
[1] "Snappp CRS:"
print(st_crs(snappp))
Coordinate Reference System:
  User input: +proj=longlat +datum=WGS84 +no_defs 
  wkt:
GEOGCRS["unknown",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
# Reprojeter au besoin (en convertisant le snappp en facteur)
snappp_sf <- st_as_sf(snappp)  

snappp_sf <- st_transform(snappp_sf, crs = st_crs(pointer3))

tm_shape(streams) +
  tm_raster(legend.show = TRUE, palette = "Blues") +
  tm_shape(snappp_sf) +
  tm_dots(col = "red", size = 0.5)
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
visual variable `col` namely 'palette' (rename to 'values') to col.scale =
tm_scale(<HERE>).
SpatRaster object downsampled to 3740 by 2675 cells.

The visual variable "col" of the layer "raster" contains a unique value. Therefore a discrete scale is applied (tm_scale_discrete).

[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
"brewer.blues"